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Brownian  bugs  and  superprocesses 

W.R.  Young 

Scripps  Institution  of  Oceanography,  University  of  California  at  San  Diego,  La  Jolla  92093-0213 


Abstract.  Advection-difFusion-reaction  (ADR.)  models  link  physical  oceanog- 
raphy and  biological  oceanography.  These  models,  which  describe  biology 
using  continuous  concentration  fields,  usually  neglect  individual-scale  fluctu- 
ations. I describe  a stochastic  individual-based  model,  called  the  Brownian 
bug  process,  which  illustrates  some  of  the  surprising  issues  associated  with 
the  neglect  of  fluctuations  by  ADR  descriptions.  The  Brownian  bug  model 
is  an  ensemble  of  random  walkers  which  suffer  birth  and  death  at  constant 
mean  rates.  (Probabilists  will  recognize  the  Brownian  bug  model  as  the 
simplest  example  of  a “superprocess.”)  Binary  division  puts  two  bugs  (parent 
and  progeny)  at  the  same  position  and  the  accumulation  of  these  small-scale 
density  fluctuations  can  produce  palpable  nonuniformities  on  large  scales. 

In  other  words,  provided  that  the  diffusion  is  not  too  strong  relative  to  the 
reproduction  rate,  a spatially  homogeneous  initial  condition  spontaneously 
develops  patches  and  voids.  The  wavenumber  signature  of  these  reproductive 
pair  correlations  is  that  the  spectrum  of  density  fluctuations  remains  white 
but  rises  linearly  with  time.  Diffusion  opposes  this  reproductive  forcing,  most 
effectively  at  large  wavenumbers,  so  that  a red  spectrum  develops. 


1.  Introduction 

Individuals  are  the  fundamental  unit  of  ecology  but 
it  is  not  feasible  to  model  ecological  processes  by  track- 
ing single  organisms.  Instead,  ecologists  often  em- 
ploy advection-diffusion-reaction  (adr)  approximations 
which  describe  biology  using  the  ‘concentration’  C(x,  t) 
(e.g.,Flierl  et  al.  1999).  The  concentration  is  defined 
via  a sample  area1  d A surrounding  any  x at  time  t: 

C(x,  t ) dA  = expected  number  of  organisms  in  dA. 

(1) 

The  dimensions  of  dA  are  such  that  the  strong  nonuni- 
formities characteristic  of  scales  comparable  to  that  of 
an  individual  are  greatly  reduced  by  averaging. 

Ecological  ADR  models  originated  with  papers  by 
Fisher  (1937),  Kolmogorov  et  al.  (1937)  and  Skellam 
(1951).  These  authors  independently  proposed 

Ct  = <yC  - A C2  + DV2C  (2) 

as  a model  for  the  growth,  saturation  and  dispersion  of 
a population.  This  chapter  of  theoretical  ecology  is  ei- 
ther a slide  down  from  this  summitt,  as  the  ADR  model 
is  applied  to  specific  cases,  or  an  ascent  as  the  ADR 


1In  this  note  I will  work  in  two-dimensions  so  that  in  (1)  I am 

using  a control  area  rather  than  a control  volume. 


approximation  is  derived  and  its  limitations  are  better 
understood.  The  1997  compilation  edited  by  Tilman 
and  Kareiva  gives  a good  overview  of  both  these  en- 
deavours. In  this  review  we  will  admire  some  scenery 
on  the  climb-up  to  the  adr  approximation  by  describ- 
ing a maximally  simplified  model — the  Brownian  bug 
process. 

Ocean  ecology  is  a natural  setting  for  ADR  descrip- 
tions because  the  spatial  domain  is  large  and  the  pop- 
ulations of  some  species  are  enormous.  Population  size 
matters  because  C in  (1)  is  an  average  or  expectation. 
The  actual  number  of  organisms  in  dA,  say  N,  is  a ran- 
dom variable  and  using  the  definition  in  (1)  N = C&A. 
The  difference  N'  = N — N is  a fluctuation.  Deter- 
ministic ADR  equations  ignore  these  fluctuations.  The 
naive  justification  for  discarding  fluctuations  is  that  N 
is  a Poisson  random  variable  (e.g.,  Durrett  and  Levin 
1994)  so  that 

N'  1 
"ST F ^ I • 

n viv 

Thus  a deterministic  description  like  (2)  may  be  useful 
for  copepods,  but  not  for  whales. 

In  oceanography,  ADR  models  are  a main  line  of  com- 
munication between  physicists  and  biologists.  Physical 
oceanographers  put  the  A into  the  ADR  equation  (2)  by 
adding  the  term  u-VC  to  the  left  hand  side.  The  ve- 
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Figure  1.  A Brownian  bug  simulation  using  A = 0.001, 
p = q — 1/4,  No  = 4000,  and  {7  = 0 (no  advection).  The 
position  of  each  bug  is  plotted  as  a point  in  the  unit  square. 
The  top  panel  shows  the  initial  condition,  the  middle  panel 
is  after  10  cycles  and  the  bottom  after  100  cycles. 


locity  u(x,  t)  is  supplied  by  some  physical  model.  (For 
example,  see  the  paper  by  Kelvin  Richards  in  these  pro- 
ceedings.) Physicists  find  this  type  of  model  appealing 
because  the  mathematical  structure  is  familiar,  and  be- 
cause it  is  easy  to  plug  biology  into  existing  ocean  mod- 
els. This  is  also  my  excuse  for  discussing  these  issues  in 
a meeting  devoted  to  oceanographic  stirring  and  mix- 
ing. 

When  I described  the  Brownian  bug  process  in  my 
presentation  at  the  2001  ’Aha  Huliko’a  meeting,  Mel 
Briscoe  pointed  out  that  the  same  model  had  been  in- 
dependently presented  at  the  1997  ' Aha  Huliko’a  meet- 
ing by  Robert  Adler  and  identified  as  a superprocess. 
Adler  approaches  this  subject  from  the  perspective  of 


a probabilist  and  the  references  in  his  article  provide 
an  entry  into  the  extensive  and  very  mathematical  lit- 
erature on  superprocesses.  A Google  search  on  “super- 
process” gives  about  250  hits,  most  of  which  are  on  the 
homepages  of  mathematicians.  Table  stakes  in  the  su- 
perprocess community  are  an  understanding  of  measure 
theory  and  branching  processes — most  oceanographers 
can’t  ante. 

Another  closely  related  class  of  problems  is  diffusion- 
limited  reactions,  e.g.,  see  the  review  by  Mattis  and 
Glasser  (1998),  which  emphasizes  the  utility  of  quan- 
tum field  theory  in  these  problems.  There  seems  to 
be  little  communication  between  physicists  interested 
in  diffusion-limited  reactions  and  probabilists  studying 
superprocesses. 

Ecologists  and  oceanographers  will  not  find  it  easy  to 
understand  the  review  articles  written  by  either  of  these 
communities.  This  review,  which  is  intended  for  my 
oceanographic  colleagues,  might  then  be  titled  “Super- 
processes for  pedestrians,  with  no  quantum  field  theory 
and  no  measure  theory.” 

2.  Brownian  bug  simulations 

The  Brownian  bug  process  is  an  ensemble  of  random 
walkers  (idealized  as  points)  milling  around  in  continu- 
ous time  and  space  and  simultaneously  reproducing  and 
dying  at  constant  mean  rates.  Reproduction  is  random 
in  time  and  occurs  via  binary  division,  so  that  parent 
and  progeny  are  at  the  same  point  for  just  an  instant 
before  their  independent  random  walks  start  to  produce 
y/Dt  separation.  Death  is  also  random  and  dead  bugs 
are  removed  from  the  ensemble.  The  addition  of  advec- 
tion to  this  scenario  ( Young  2001)  results  in  the  simplest 
model  of  a planktonic  species  reproducing  and  dying  in 
a turbulent  fluid. 

To  simulate  the  continuous  process  described  above 
on  a computer  one  must  discretize.  There  are  sev- 
eral ways  of  doing  this,  and  I describe  in  the  appendix 
an  individual-based  Monte  Carlo  procedure  which  pro- 
duces Figure  1.  The  top  panel  shows  the  homogeneous 
initial  condition  and  the  middle  and  bottom  panels  the 
results  after  10  and  100  cycles  of  birth,  death  and  dif- 
fusion. Patches  of  bugs  spontaneously  form  out  of  the 
spatially  uniform  initial  conditions. 

Failure  of  deterministic  ADR.  The  ADR  ap- 
proximation of  the  simulation  in  Figure  1 is  simply 

Ct  = DV2C+(\-p)C,  (4) 

where  D is  the  diffusivity,  A is  birth  rate,  and  p the 
death  rate.  (For  the  relation  between  D,  A,  and  p and 
the  parameters  in  the  simulation  see  the  appendix.)  But 
in  Figure  1 birth  and  death  are  equiprobable,  X = p, 
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and  consequently  (4)  collapses  to  the  diffusion  equa- 
tion. Thus,  with  the  initial  condition  C(x,  0)  = Nq/L2, 
the  solution  of  the  ADR  description  is  C(x,  t)  = Nq/L2. 
Uniform  concentration  is  not  a good  characterization  of 
Figure  1 and  we  conclude  that  the  ADR  approximation 
fails  to  describe  the  Brownian  bug  process. 

A stochastic  ADR  equation.  Although  deter- 
minstic  ADR  fails,  there  is  a stochastic  partial  differen- 
tial equation  model  which  does  capture  the  fluctuations 
exhibited  by  the  individual-based  model.  This  stochas- 
tic equation  is  something  like 

Ct  = DV2C  + VCw(x,t)  (5) 

where  w is  a white  Gaussian  noise.  The  yfC  in  (5) 
occurs  because  in  the  control  volume  d the  stochastic 
imbalance  between  births  and  deaths  is  proportional  to 
y/C.  I qualified  the  introduction  of  (5)  with  “something 
like”  because  mathematicians  (e.g.,  Adler  1997)  tell  us 
that  (5)  makes  sense  only  in  one  dimension,  d — 1,  but 
not  when  d>  2:  in  that  case  the  concentration  C(x,  t) 
cannot  be  defined.  I do  not  fully  understand  this  point. 
In  any  event,  Adler  proceeds  to  introduce  an  integrated 
version  of  (5)  which  he  claims  is  sensible  even  in  d > 2. 
We  press  on... 

Density  fluctuations.  The  spatial  clumping  in 
Figure  1 occurs  because  birth  is  always  next  to  a living 
bug,  while  death  can  occur  anywhere.  Alternatively,  in 
some  regions  because  of  bad  luck  there  is  local  extinc- 
tion of  the  bug  population.  Diffusion  from  the  surviving 
centers  of  population  is  not  fast  enough  to  recolonize 
the  voids.  This  suggests  that  to  produce  patches  the 
diffusion  D must  be  weak  relative  to  the  growth  rate 
A = p.  We  can  identify  the  nondimensional  parameter 
that  controls  this  transition  by  noting  that  the  initial 
separation  distance  between  bugs  (the  “specific”  length) 
is  Is  = L/y/No . On  the  other  hand,  using  the  diffusivity 
D and  the  growth  rate  A,  we  can  define  a “reproduc- 
tive” length  iT  — y/D/X,  which  measures  the  separa- 
tion that  develops  between  parent  and  progeny  before 
there  is  another  reproductive  event.  We  expect  that 
patchiness  will  be  strong  if  the  nondimensional  number 


Figure  2.  Results  of  7 simulations  each  with  No  = 218 
bugs  at  t = 0.  The  unit  square  is  divided  into  a lattice 
of  32  x 32  boxes  (i.e.,  p = 5)  and  <r(5,  t)  defined  in  (10)  is 
computed.  The  parameters  in  section  5 are  p = q = 1/4 
and  the  step-length  of  the  random  walk,  A,  is  varied  as 
indicated.  When  the  diffusion  is  strong  <r  remains  close 
to  its  initial  Poisson  value,  er(5,0)  « 1/16.  The  dotted 
curve  is  ctgw  = y/{l  + 2 pt)/np,  which  is  computed  using 
the  Galton- Watson  solution  in  (18)  and  (19)  below. 


There  are  22p  boxes  so  that  the  expected  number  of 
bugs  per  box  is 

np  & 2~2pN0 . (8) 

The  occupation  number  of  box  b is  the  fluctuating  quan- 
tity 


rib(t)  = the  number  of  bugs  in  box  b at  t,  (9) 

where  b = 1, . . . , 22p  indexes  the  boxes.  To  measure  the 
size  of  fluctuations  in  rib  relative  to  the  expected  value 
rip  we  use  the  ratio  of  the  RMS  value  of  n(b,  t ) to  np: 


a(p,t) 


22p 


N 


]T[u(m) 


(10) 


is  large.  If  the  ratio  ls/£r  is  small  then  density  fluc- 
tuations should  be  small  because  the  diffusion  will  be 
effective  at  mixing  a family  into  the  background  of  un- 
related, and  therefore  uncorrelated,  bugs. 

To  illustrate  these  considerations,  the  LxL  square  is 
divided,  like  a chessboard,  into  smaller  squares  (“boxes”), 
each  of  size 

£p  = 2~PL.  (7) 


If  <r(p,  f)Cl  then  n(b,  t)  is  fluctuating  weakly  around 
the  expected  value  np. 

Figure  2 shows  a(p,  t)  for  a set  of  simulations  with 
various  diffusivities  and  p = 5.  The  a(p,  0)  is  small  be- 
cause the  box  size  tp  = 1/32  is  much  greater  than  the 
specific  length,  ls  = 1/512.  If  the  diffusivity  is  strong 
enough  (for  example,  A > 2-8  in  Figure  2)  then  a(jp,  t) 
experiences  mild  growth  but  cr(p,  oo)  -C  1.  In  this  case 
the  diffusivity  is  arresting  reproductively  driven  fluctu- 
ations. On  the  other  hand,  if  the  diffusion  is  small  then 
the  growth  of  a proceeds  till  a(p,oo)  > 0(1),  which 
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indicates  large  fluctuations  in  the  density  defined  using 
the  box-scale  lp.  In  particular,  if  there  is  no  diffusion 
(A  = 0)  then  one  can  show  using  the  Galton- Watson 
solution  in  (18)  and  (19)  that  a(p,t)  « ^/(l  + 2 pt)/np. 
Notice  that  by  decreasing  the  diffusivity  we  are  increas- 
ing the  nondimensional  parameter  in  (6). 

It  would  be  useful  to  be  able  to  give  a scaling  law  for 
the  saturated  value,  a(p,  oo),  in  terms  of  D,  p,  lp  and 
the  other  external  parameters.  I have  not  been  able  to 
achieve  this. 

3.  Galton- Watsonology 

The  Brownian  bug  model  combines  two  classical 
stochastic  processes:  Brownian  motion  and  a Galton- 
Watson  branching  process.  Brownian  motion  is  familiar 
to  most  oceanographers,  but  the  Galton- Watson  pro- 
cess is  not.  This  section  summarizes  some  basic  prop- 
erties of  Galton- Watson  process. 

Suppose  that  we  switch  off  the  motion  of  the  bugs, 
i.e.,  we  set  the  step  length  of  the  random  walk  to  zero. 
In  this  case  each  bug  founds  a family  which  sits  at  the 
initial  spatial  location.  Newly  born  bugs  pile  up  and 
build  a tower  on  top  of  their  ancestors.  The  towers  go 
up  and  down  as  bugs  are  born  and  die.  Towers  can  dis- 
appear forever  (extinction)  if  a run  of  bad  luck  wipes 
out  the  family.  The  statistics  of  these  families,  particu- 
larly the  probability  of  extinction,  were  studied  by  Gal- 
ton and  Watson  in  the  nineteenth  century  (see  Harris, 
1963,  and  chapter  XVII  of  Feller , 1968).  In  this  review 
we  closely  follow  Feller’s  continuous-time  formulation  of 
the  Galton- Watson  process. 

The  fundamantal  variables  of  the  Galton- Watson 
process  are  the  probabilities 

Pk(t ) = the  probability  that  there  are  k bugs  at  t. 

We  study  the  proliferation  of  a family  descended  from 
a single  urbug  using  the  initial  condition  p\ (0)  = 1 and 
Pk(0)  = 0 if  k 7^  1.  The  probability  that  a family  is 
extinct  at  t is  po(t).  Given  pfc(t),  one  calculates  the 
average  of  a function  of  k according  to  the  definition 

OO 

/(M)  = ^/(M)Pfc' (t).  (11) 

fc= o 

Normalization  of  probabilities  requires  1 = 1 and  the 
mean  population  of  a family  is  k(t). 

The  family  becomes  extinct  when  a family  with  one 
surviving  member  suffers  a death;  accordingly  the  prob- 
ability of  extinction  po  increases  like 

(12) 


The  evolution  of  the  other  probabilities  with  k > 1 is 
obtained  from  the  hierarchy  of  differential  equations 

pk  = A (k  - l)pjt— i - (A  + p)kpk  + p(fc  + l)pfc+i.  (13) 

The  term  — (A  + p)kpk  on  the  right  hand  side  of  (13) 
appears  because  a family  of  k bugs  might  make  a tran- 
sition to  k + 1 bugs  if  any  of  its  k members  repro- 
duce (probability  per  time  A k)  or  make  a transition  to 
k — 1 bugs  if  any  of  the  k members  die  (probability  per 
time  pk).  Analogous  reasoning  is  used  to  deduce  the 
other  terms  containing  k — 1 and  k + 1 on  the  right 
hand  side  of  (12). 

Mean  and  variance.  From  (12)  and  (13)  one  can 
show  using  the  definition  (11)  that  1 = 1 and 

ft=(\-p)k.  (14) 

This  result  establishes  the  connection  between  the  Gal- 
ton-Watson  process  and  the  ADR  equation  (3):  both 
k and  C are  expectations  defined  by  ensemble  averages 
and  A - p is  the  growth  rate  of  the  expected  population. 
In  fact,  if  we  start  with  N0  bugs,  each  of  which  is  a 
family  with  one  member,  then  because  the  evolution  is 
independent  f CAA  = N0k. 

The  statistic  k contains  no  information  concerning 
the  fluctuations  in  population.  That  is,  at  times  t > 0 
some  families  are  extinct,  while  others  are  very  large.  A 
coarse  description  of  these  fluctuations  is  provided  by 
the  variance  k2  -k2,  or  better  by  the  ratio 


v = Vk2  — k2/k , (15) 

which,  like  (3),  compares  the  fluctuations  to  the  mean. 
Assiduous  summation  of  (12)  and  (13)  shows  that 

dk2 

— = 2(A-p)fc2  + (A  + p)fc.  (16) 

at 

If  the  death  rate  balances  the  birth  rate,  so  that  A = p 
and  k = 1,  then  the  solution  of  (16)  is 

k2  = 1 + 2pf , (A  = p) . (17) 

In  this  case  v = %/2pt,  indicating  that  the  system  is 
dominated  by  fluctuations  a sf-+oo. 

The  special  case  A = p in  (17)  is  the  critical  point 
of  the  branching  process.  If  7 = A — p>0  then 
k ~ exp(7f)  — ^ 00  and  v ~ exp(-7t/2)  ->  0.  On 
the  other  hand,  if  7 < 0 then  v ~ exp(— 7f/2)  — > 00. 
The  conclusion  is  that  the  system  is  dominated  by  fluc- 
tuations if  the  growth  rate  7 is  less  than  or  equal  to 
zero. 


po  = m • 
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A solution  with  A = p.  The  variance  is  only  a 
rough  guide  to  behaviour  of  fluctuations.  A much  more 
informative  result  is  provided  by  an  explicit  solution  of 
(12)  and  (13)  which  is  obtained  using  the  generating 
function  method  ( Feller  1968): 

= (18) 


and  if  k > 1: 


Pfc(t)  = 


(1  + pt)k+1  ' 


(19) 


At  large  times,  pt  2>  1,  see  that  most  families  are  ex- 
tinct, po(t)  — > 1.  That  is,  at  large  times  the  most  prob- 
able number  of  members  in  any  given  family  is  zero  — 
in  fact  extinction  is  certain.  On  the  other  hand  k = 1 
so  that  the  expected  number  of  members  in  a family 
(counting  extinct  families  as  having  zero  members)  is 
always  constant. 


To  reconcile  the  extinction  of  most  families  with  a 
constant  expected  population,  we  make  a distinction 
between  surviving  families  (meaning  a family  with  at 
least  one  member)  and  extinct  families.  Suppose  that 
we  start  with  a large  number  N0  » 1 of  individuals, 
each  of  which  founds  a family.  At  t the  number  of  sur- 
viving families  is  a random  variable,  F(t).  The  expected 
number  of  surviving  families  is  F = (1  —po)N0,  or  using 
(18) 


N0 

1 + pt' 


(20) 


Let  the  random  variable  n be  the  number  of  individuals 
in  a surviving  family.  Prom  (18)  and  (19),  the  proba- 
bility distribution  of  n,  conditioned  on  the  survival  of  a 
family,  is 


(Ft)n 

(1  + pt)n  ’ 


(n  = 1, 2 • • •) . 


and  consequently 


(21) 


4.  Spectral  dynamics 

We  began  section  3 by  imagining  that  the  bugs  were 
immobile  so  that  the  surviving  families  formed  towers 
on  top  of  the  urbug.  We  then  ignored  the  spatial  aspects 
of  the  model,  and  focussed  simply  on  the  statistics  of  the 
bug  population.  Now  let  us  return  to  the  spatial  domain 
and  study  the  spectrum  of  concentration  fluctuations. 

For  the  moment,  we  continue  to  assume  that  the  bugs 
are  immobile  ( D = 0).  The  discussion  surrounding 

(20)  through  (22)  shows  that  the  height  of  a typical 
family  tower  grows  like  1 + pt  and  the  average  sep- 
aration between  towers  grows  like  £s/y/l~+pt,  where 
£s  = L/y/Wo-  We  can  write  the  microscopic  density  of 
bugs  in  a single  realization  of  the  process  as 

F 

p(x,  t)  = Y n/<*(x  - x/)  > (23) 

/= i 

where  F(t)  is  the  number  of  surviving  families,  n/  is  the 
number  of  members  in  family  / and  X/  is  the  location 
of  the  family.  Notice  that  the  concentration  is 

C = (p)  = N0/L2 , (24) 

where  the  ()  denotes  an  ensemble  average  over  the  ran- 
dom locations  2 X/  and  an  average  over  the  n/’s  using 

(21) .  The  average  C = (p)  is  trivial  — just  like  k in  the 
previous  section.  To  learn  something  about  the  fluctu- 
ations we  must  study  quadratic  quantities  such  as  the 
spectrum. 

Continuing  with  the  D = 0 assumption,  we  now 
make  a standard  calculation  of  the  spectrum  by  ensem- 
ble averaging  the  density  in  (23).  The  Fourier  series 
representation  of  (23)  is 

p{x,  t)  = L~2  Y Pk  exp  (ik-  x)  , (25) 

k 


where 


F 

Pk  = ’Y,  nf  exp  (-ik-x/)  . 

/= i 


(26) 


n(f)  = pt  + 1 , n2  = (2/it  + 1 )(pt  + 1) . (22) 

Thus,  the  number  of  surviving  families  decreases  like 
t-1,  while  the  membership  of  a surviving  family  grows 
like  t 1 to  keep  the  total  population  fixed. 

These  remarks  illustrate  Galton  and  Watson’s  con- 
clusion that  the  “observed  diminution  of  surnames  among 
families  whose  history  we  can  trace,  is  not  a sign  of  their 
diminished  fertility.”  Instead  the  extinction  of  most 
families,  and  the  apparent  fecundity  of  a few  exceptions, 
is  a result  of  “the  ordinary  law  of  chances.” 


The  spectrum  is  related  to 

(pkPk) - (Yn* +2Y,Ysnfnf'^k{x,'~Xf)) ’ (27) 

/= i f>f 

where  the  ()  is  an  ensemble  average  over  both  the  ran- 
dom and  uncorrelated  positions  x/  of  the  surviving  fam- 
ilies and  over  the  random  variable  nj.  Because  there  are 


2The  family  locations  Xy  are  uniformly  distributed  in  the  do- 
main. Thus  (S(x  — Xf))  = L~2  and  if  k ^ 0 then  (exp(ik-xy))  = 
0. 
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Figure  3.  Comparison  of  the  analytic  prediction  (31)  (the 
smooth  solid  curve)  with  the  results  of  simulation.  We  show 
estimated  spectra  (the  three  dotted  curves)  from  the  simula- 
tion at  t = 1,  t = 10  and  t = 100;  the  comparison  with  (31) 
is  made  only  at  t = 100.  The  simulation  used  No  = 105  bugs 
with  A = 2-u  and  V = q~  7/64.  The  spectra  are  estiamted 
by  azimuthally  averaging  in  the  wavenumber  plane. 

no  correlations  between  the  positions  of  the  families  the 
final  term,  involving  the  double  sum  with  / /'  van- 

ishes. Using  F from  (20),  and  n2  from  (22),  shows  that 
the  remaining  term  on  the  right  hand  side  of  (27)  is 

(PkPu)  = Fn2  = JVb(l  + 2 pa)  - (28) 

This  simple  calculation  indicates  that  the  spatial  sig- 
nal of  a diffusionless  Galton- Watson  process  is  a white 
spectrum  which  rises  at  a constant  rate. 

How  is  the  conclusion  above  affected  by  diffusion? 
With  D 0 the  members  of  a family  move  as  indepen- 
dent random  walkers  and  the  family  towers  blur  into 
spreading  clouds  of  bugs.  We  could  write  down  a micro- 
scopic density,  analogous  to  (23),  in  which  the  position 
of  each  individual  bug  is  represented  with  a 5-function. 
But  then  the  calculation  founders  at  (27)  because  the  {) 
no  longer  kills  the  off-diagonal  terms:  there  are  correla- 
tions between  the  positions  of  bugs  which  belong  to  the 
same  family.  Instead  we  can  argue  heuristically  that 
the  evolution  of  the  spectrum,  defined  as 

S(k,t)  = N^(PkPl),  (29) 

is  given  by 

St  = 2p-  2 Dk2S , S(k,  0)  = 1 , (30) 

where  k = Vk  k.  The  first  term  on  the  right  hand 
side  of  (30)  is  a source  of  variance  which  produces  the 
uniform  rise  of  the  spectral  level  in  (28).  The  the  final 


term  is  the  diffusive  sink  of  variance  which  ultimately 
balances  the  source. 

The  solution  of  (30) 

S{Kt)  = DV  + (31) 

exhibits  a k~~ 2 range  at  high  wavenumbers.  (If  we 
isotropize  by  using  2nkS(k ) then  this  is  a spec- 
trum.) Figure  3 shows  good  comparison  of  (31)  with  a 
spectrum  estimated  from  a simulation. 

5.  Conclusions 

In  this  note  we  have  studied  a model  which  focuses 
attention  on  the  role  of  fluctuations  in  the  spatial  de- 
velopment of  populations.  It  is  remarkable  that  Euler, 
Navier  and  Stokes  formulated  the  basic  equations  of 
fluid  mechanics  before  the  reality  of  atoms  was  estab- 
lished. This  history  indicates  that  for  most  purposes 
atomic  fluctuations  are  unimportant  in  fluid  mechanics. 
The  Brownian  bug  model  suggests  that  the  problem  of 
biological  fluctuations  is  more  severe,  and  that  events 
on  microscopic  scales  can  impact  macroscopic  observa- 
tions. Indeed,  the  source  term  2fi  on  the  right  hand 
side  of  (30)  is  due  to  reproduction  and  has  no  analog  in 
fluid  mechanics. 

Appendix:  A recipe  for  simulations 

The  initial  condition,  t = 0,  is  prepared  with  a large 
number  No  of  bugs,  strewn  randomly  into  an  L x L 
square  domain.  The  initial  condition  is  a Poisson  point 
process  in  which  there  are  no  correlations  between  the 
positions  of  the  bugs. 

The  simulation  is  advanced  through  time  in  incre- 
ments of  a “cycle”  of  duration  r.  Each  cycle  consists 
of  two  steps:  (a)  random  birth  and  death;  (b)  Brow- 
nian motion.  In  step  (a)  each  bug  reproduces  by  bi- 
nary fission  (probability  p)  or  dies  (probability  q)  or 
remains  unchanged  (probability  1 — p — q).  When  a 
lucky  bug  divides  the  offspring  is  placed  on  top  of  the 
parent.  In  step  (b),  bug  k is  displaced  to  a new  position 
x*;(f  + t)  = Xfc(t)  + 5xjt(t).  The  components  of  Sxk  are 
independent  and  identically  distributed  Gaussian  ran- 
dom variables,  each  with  RMS  value  A (i.e.,  Brownian 
motion  with  the  diffusivity  D = A2/2r).  These  inde- 
pendent displacements  separate  the  coincident  parent- 
progeny  pairs  created  by  step  (a). 

In  principle,  one  can  approach  the  continuous  limit 
by  taking  r ->  0 while  holding  the  parameters  A = p/r, 
fi  - q/r,  and  D = A2/2t  fixed. 

We  can  also  use  po{t)  in  (18)  to  calculate  the  prob- 
ability that  an  initial  population  of  No  bugs  becomes 
extinct  at  time  t : each  bug  produces  a family  which 
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evolves  independently  of  the  other  families,  so  that  the 
probability  that  all  Nq  families  are  extinct  at  t is 

Pext  = Po°  ~ exp (-Text/t)  , (Al) 

where  the  extinction  time-scale  is  Text  = Nq/h-  (The 
exponential  approximation  in  (Al)  assumes  that  No  is 
large  and  t/Text  is  of  order  unity.)  In  simulations  we 
avoid  extinction  by  making  No  much  larger  than  the 
number  of  generations.  The  suppression  of  extinction 
is  a consequence  of  approaching  the  “thermodynamic 
limit”  in  which  JV0  -»  oo  with  Nq/L2  fixed. 
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